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ABSTRACT 

We report the detection of three transiting planets around a Sunlike star, which we designate Kepler- 
18. The transit signals were detected in photometric data from the Kepler satellite, and were confirmed 
to arise from planets using a combination of large transit-timing variations, radial- velocity variations, 
Warm-Spitzer observations, and statistical analysis of false-positive probabilities. The Kepler-18 star 
has a mass of O.97M0, radius I.IRq, effective temperature 5345K, and iron abundance [Fe/H]= +0.19. 
The planets have orbital periods of approximately 3.5, 7.6 and 14.9 days. The innermost planet "b" 
is a "super-Earth" with mass 6.9 ± 3.4M®, radius 2.00 ± 0.10-R®, and mean density 4.9 ± 2.4g cm 3 . 
The two outer planets "c" and "d" are both low-density Neptune-mass planets. Kepler-18c has a 
mass of 17.3 ± 1.9M®, radius 5.49 ± 0.26i? e , and mean density 0.59 ± 0.07g cm 3 , while Kepler-18d 
has a mass of 16.4 ± 1.4M®, radius 6.98 ± 0.33-R®, and mean density 0.27 ± 0.03 g cm 3 . Kepler-18c 
and Kepler-18d have orbital periods near a 2:1 mean-motion resonance, leading to large and readily 
detected transit timing variations. 

Subject headings: planetary systems — stars: individual (Kepler-18, KIC 8644288, 2MASS J19521906 
+4444467) — techniques: photometric — techniques: spectroscopic 
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1. INTRODUCTION 

Kepler is a NASA Mission designed to detect the tran- 
sits of exoplanets across the disks of their stars. The 
ultimate mission goal is to detect the transits of po- 
tentially habitable Earth-size planets. To achieve this 
goal requires a telescope in a very stable space environ- 
ment with a large (0.95-meter) effective aperture moni- 
toring the brightness of about 150,000 stars simultane- 
ously and continuously for over three years. The Ke- 
ple r Mission design and perform ance are summarize d 
by iBorucki et all (|2010bl ) and by iKoch et al.1 (|2010bt l. 
and a discussion of the commiss i oning a n d first quarter 
data ar e given by IBorucki et al.l (|2011al ). IBorucki et al.l 
(|2011bh reported 1235 planet candidates that were dis- 
covered during the fi rst four months of the Mission. 
iBatalha et al.l (|2010bD discuss the selection and charac- 
teristics of the Kepler target stars. 

The first 5 planets discovered by the Kepler mis- 
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sion (Kepler 4-8) were reported in January 2010 
(IBorucki et al.ll2010a| iKoch et al.lfeOlOal [Dunham et al.l 



120101: iLatham et al.M2010t Uenkins et alil2010al) . Ke- 
pler has detec t ed an a bundance of multi-planet systems. 
IBorucki et al.l (|2011bl ) reported a total of 170 candi- 
date multi-planet systems among the 997 planet can- 
didate host stars f rom the February 2011 data release. 
ISteffen et al.l (|2010f ) presented five o f" these systems in de- 
tail. The Kepler-9b and c system ([Holman et al.l 120101) 
was the first transiting multi-planet sys tem confirmed 
by tr ansit timing variations. Kepler-10 (Batalha ct al. 
Honj) was the first rocky p lanet found by Kepler. Kepler- 
11 ( Lissauer et al.ll2011al ) is a transiting system of six 
planets. 

Not all Kepler planets can be directly confirmed by 
supporting reflex radial velocity measurements of the 
parent star, or by detection and modeling of transit tim- 
ing variations. Instead, some planets must be validated 
by analyzing all possible possible astrophysical false- 
positive scenarios and comparing their a priori likelihood 
to that of a planet. The Kepler project has been able to 
utilize the BLENDER technique developed by iTorres et al.l 
(l2004fl to validate a t hird planet in the Kepler-9 sys- 
tem jTorres et alll201lD . a second planet in the Kepler-10 
(|Fressin et a,l.ll2011|) and the out er planet in the Kepler- 
11 system ([Lissauer et al.ll2011a| ). 

Here we present the Kepler-18 system, containing two 
Neptune-mass transiting planets near a 2:1 mean mo- 
tion resonance which show significant gravitational inter- 
actions which are observed via measurements of transit 
timing variations (TTVs) , as well as a small, inner super- 
Earth size transiting planet. This system is remarkably 
similar to the Kepler-9 system in its overall architecture. 

2. Kepler PHOTOMETRY 

The Kepler spacecraft carries a photometer with a 
wide-field (~ 115 deg 2 ) Schmidt camera of 0.95-m effec- 
tive aperture. The spacecraft was launched in March 
2009, and is now in an Earth-trailing heliocentric or- 
bit which allows nearly continuous photo metric coverage 
of its field-of- view in Cygnus and Lyra. ICaldwell et al.1 
(2010) discuss the early instrumental performance of the 
Kepler photometer system. The primary data for de- 
tection of transiting planets are the Long Cadence (LC) 
"Pre-search Data Conditioned" (PDC) time series data, 
in which 270 consecutive CCD readouts are binned, giv- 
ing an effective sam pling interval of 29.4244 minutes 
([Jenkins et al.l l2010bl ). A small selected subset of Ke- 
pler targets is sampled at the Short Cadence (SC) rate 
of 9 c onsecutive reads for a sampling interval of 58.85 sec- 
onds ([Gilliland et aLll2010h . Thus, one LC sample is the 
sum of 30 SC samples. The data from the spacecraft are 
processe d through the Kepler Science Operations Center 
pipeline ([Jenkins et al.ll2010cl ) to perform standard CCD 
processing and to remove instrumental artifacts. The 
LC PDC time series data are searched for possible plan- 
etary transits using a wav elet-based adaptive matched 
filter ( J enkins et al.l l2010dl ). Possible planetary transit 
events with amplitude greater than 7.1a are flagged and 
are then subjecte d to intensive valida t ion efforts usin g 
the Kepler data ([Batalha et al.ll2010at l"Wu et al.ll2010h . 
Objects that pass this level of vetting are designated as 
a "Kepler Object of Interest" (KOI) and are sent to the 
Follow-up Observing Program (FOP) for further study. 



2.1. Light Curves and Data Validation 

One of the objects identified with possible transiting 
planets is the Kp 13.549 magnitude (where Kp is the 
magnitude in the Kepler passband) star KIC 8644288 
(2MASS J19521906+4444467, K00137). After a possible 
transiting planet has been detected, the Kepler data are 
subjected to a set of statistical tests to search for possible 
astrophysical false-positive origin of the observed signal. 
These data validation tests f or the first five Kepler planet 
discoveries are described by iBatalha et al.l (|2010aD . Ad- 
ditional tests, including meas urement of the image cen- 
troid motion during a transit (|Wu et al.ll2010D all gave a 
high probability that the signals seen were real. The ap- 
plication of these techniques to Kepler- 10b is described 
in detail bv IBatalha et all ([201 If ). 

Two separate transiting objects were immediately ob- 
vious in the LC data. K00137.01 has a transit ephemeris 
of T [BJD] = (2455167.0883 ± 0.0023) + N * (7.64159 ± 
0.00003) days and a transit depth of 2287 ± 9 ppm. (All 
transit times and ephemerides in this paper are based 
on UTC.) K00137.02 has an ephemeris of T [BJD] = 
(2455169.1776 ± 0.0013) + N * (14.85888 ± 0.00004) days 
and a transit depth of 3265 ± 12 ppm. It was noted that 
the orbits of these two objects were very near a period 
ratio of 2:1. After filtering these transits of K00137.01 
and K00137.02 from the lightcurve, we searched again for 
transiting objects, and found a third planet candidate in 
the system, K00137.03, which has a shorter orbital pe- 
riod than the other two transiting planets. K00137.03 
has the ephemeris T [BJD] = (2454966.5068 ± 0.0021) + 
N * (3.504725 ± 0.000028) days, and a transit depth of 
only 254 ± 8 ppm. The light curve for K00137 is shown 
in Figure [TJ The upper panel shows the raw uncorrected 
"Photometric Analysis" (PA) lightcurves that come out 
of the data processing pipeline, and the lower panel 
shows the corrected PDC lightcurves. The two tran- 
sit events that look significantly deeper than the others, 
near BJD 2454976.0 and BJD 2455243.5 are simultaneous 
transits of K00137.01 and K00137.02. 

The folded lightcurves for each of the three candidate 
planets are shown in Figure [2j This Figure shows the 
lightcurves folded on the ephemeris given above for each 
KOI. The significant width of the ingress and egress for 
Kepler-18c and Kepler-18d are a result of the transit time 
variations discussed in Section [6l 

3. FOLLOW-UP OBSERVATIONS 

After the possible transiting planets were found, and 
the KOIs passed the Data Validation tests for false posi- 
tive signals, K00137 was sent on to the Kepler Follow-up 
Observing Program (FOP) for ground-based telescopic 
observations designed either to find any additional indi- 
cation that these KOIs might be an astrophysical false- 
positive signal, or to verify the planetary nature of the 
transit events. 

3.1. Reconnaissance Spectroscopy 

The first FOP step is to obtain high spectral resolu- 
tion, low S/N spectroscopy in order to verify the Ke- 
pler Input Catalog stellar classification, and to search 
for evidence of stellar multiplicity in the spectrum. Spec- 
tra from the Hamilton echelle spectrograph on the Lick 
Observatory 3-m Shane Telescope were obtained on the 
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Fig. 1. — Kepler lightcurves for K00137. The upper panel shows the normalized raw Photometric Analysis (PA) lightcurves for quarters 
through 6. The various discontinuities are due to effects such as spacecraft safe-mode events or loss of fine pointing. Each quarter put 
the star on a different detector, which accounts for the change in overall sensitivity. The long-term drifts in each quarter are temperature 
related. The lower panel shows the normalized corrected Pre-search Data Conditioned (PDC) lightcurve. Most of the spacecraft-related 
variability of the PA lightcurve has been removed. 



nights of 8 and 9 August 2009, and on 1 September 2009 
UT. These spectra showed no convincing evidence for 
radial velocity variability at the 0.5 km s -1 level, and 
no hints of any contaminating spectra. These spectra 
were cross-correlated again st a library of s ynthe tic stel- 
lar spectra as described by iBatalha et al.l ()201lD . in or- 
der to derive basic stellar parameters to compare with 
the Kepler Input Catalog (KIC) values. These spec- 
tra yielded T eff = 5250 ± 125K, logg = 4.0 ± 0.25, and 
Vsini = 2±2kms~ 1 . The height of the cross-correlation 
peaks ranged from 0.82 to 0.89, indicating a very good 
match with the library spectra. 

3.2. High Spatial- Resolution Imaging 

High resolution imaging of the surroundings of a KOI 
is an important step to identify possible sources of false- 
positive signals. We need to ensure that the detected 
transit signal is indeed originating on the selected target 
star, and not on a background star that was unresolved 
in the original KIC imaging of the Kepler field. 



A seeing limited image was obtained at Lick Obser- 
vatory's 1-m Nickel telescope using the Direct Imaging 
Camera. This image was a single one-minute exposure 
taken in the I-band, with seeing of approximately 1.5". 
The only nearby star is a faint object (approximately 
5 magnitudes fainter than K00137) about 5.5" to the 
north. A J-band image of the 1' x 1' field of view around 
K00137was taken as part of a complete J-band survey 
of the Kepler field of view using the wide field camera 
(WFCAM) on the United Kingdom Infrared Telescope 
(UKIRT). These images have a typical spatial resolution 
of 0.8-0.9" and an image depth of J = 19.6 (Vega sys- 
tem). This image also shows only a faint source about 
5.5" to the north. 

We perform speckle observations at the WIYN 3.5-m 
telescope, usi ng the Differential Speckle Survey I nstru- 
ment (DSSI) porch et alJl2011tlHowell"et^[2Ql . The 
DSSI provides simultaneous observations in two filters 
by employing a dichroic beam splitter and two identical 
EMCCDs as the imagers. We generally observe simulta- 
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Fig. 2. — Folded Lightcurves for K00137. The top row is Kepler- 
18b (K00137.03), the middle row is Kepler-18c (K00137.01) and 
the bottom row is Kepler-18d (K00137.02). The lightcurves are 
folded on the mean period listed in Table [T] 




Fig. 3. — High spatial resolution A/O image of K00137 taken 
with the Palomar Observatory 5m adaptive optics system. The 
nearest contaminating source is about 5.5" to the north. 

neously in V and R bandpasses, where V has a central 
wavelength of 5620A, and R has a central wavelength of 
6920A, and each filter has a FWHM of 400A. The speckle 
observations of K00137 were obtained on 22 June 2010 
UT and consisted of three sets of 1000, 80 msec individ- 
ual speckle images. Along with a nearly identical V-band 
reconstructed image, the speckle results reveal no com- 
panion star near K00137 within the annulus from 0.05" 
to 1.8" to a limit (5<r) of 4.2 magnitudes fainter in both 
V and R than the Kp 13.549 target star. 

Near-infrared adaptive optics imaging of K00137 was 
obtained on the night of 08 September 2009 UT with 
the Palomar Hal e 5-m telescope and the PHARO near- 
infrared camera (|Havward et al.ll200lD behin d the Palo- 
mar adaptive optics system (|Trov et al.ll2000f ) . PHARO, 
a 1024 x 20124 HgCdTe infrared array, was utilized in 
25.1 mas/pixel mode yielding a field of view of 25". Ob- 
servations were performed in J filter. The data were 



TABLE 1 

Palomar AO source sensitivity as a 

FUNCTION OF DISTANCE FROM THE 
PRIMARY TARGET AT J. 



Distance 


Distance 


AJ 


J 


(FWHM) 


(") 


(mag) 


(mag) 


1 


0.075 


1.5 


13.7 


2 


0.150 


2.5 


14.7 


3 


0.225 


3.5 


15.7 


1 


0.300 


4.0 


16.2 


5 


0.375 


4.5 


16.7 


6 


0.450 


5.0 


17.2 


7 


0.525 


6.0 


18.2 


8 


0.600 


7.0 


19.2 


9 


0.675 


7.5 


19.7 


40 


3.000 


8.5 


20.7 



collected in a standard 5-point quincunx dither pattern 
of 5" steps interlaced with an off-source (60" East) sky 
dither pattern. Data were taken with integration times 
per frame of 30 sec with a total on-source integration time 
of 7.5 minutes. The adaptive optics system guided on the 
primary target itself and produced a central core width 
of FWHM = 0.075". Figure [3] shows this Palomar im- 
age of K00137. There are two additional sources within 
15" of the primary target. The first source, located 5.6" 
to the north of K00137, is 4 magnitudes fainter at J; the 
second, located 15.0" to the southeast, is 5.8 magnitudes 
fainter. To produce the observed transit depths in the 
blended Kepler photometry, the Sj = 4 mag star would 
need to have eclipses that are 0.09 and 0.5 magnitudes 
deep, and the SJ = 5.8 mag star would need to have 
eclipses that are 0.15 and 1.3 magnitudes deep. Such 
deep eclipses from stars separated by more than 4" (> 1 
Kepler pixel) would easily be detected in the centroid 
motion analysis, but no centroid motion was detected 
between in and out of transit, indicating that these two 
stars were not responsible for the observed events. Six 
additional sources were detected at J within 15" of the 
primary target. But all of these sources are S J > 8 mag- 
nitudes fainter than K00137 and could not produce the 
observed transit events (i.e., even if the stars dimmed by 
100%, the resulting transit in the blended photometry 
would not be deep enough to match the observed transit 
depths). 

Source detection completeness was evaluated by ran- 
domly inserted fake sources of various magnitudes in 
steps of 0.5 mag and at varying distances in steps of 
1.0 FWHM from the primary target. Identification of 
sources was performed both automatically with the IDL 
version of DAOPhot and by eye. Magnitude detection 
limits were set when a source was not detected by the 
automated FIND routine or was not detected by eye. 
Within a distance of 1 — 2 FWHM, the automated find- 
ing routine often failed even though the eye could discern 
two sources, but beyond that distance the two methods 
agreed well. A summary of the detection efficiency as a 
function of distance from the primary star is given Ta- 
ble ED 

A major source of false positive planet indication in 
the Kepler data is background eclipsing binary (BGEB) 
stars within the photometric aperture of Kepler-18, 
which, when diluted by Kepler-18 itself, can produce a 
planetary-size transit signal. We perform a direct mea- 
surement of the source location via difference images. 
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TABLE 2 
Mean Pixel Response Function Fit 
Source Offsets 



Planet 


Distance (arcsec) 


a 


Kepler-18b 


0.536 ± 0.245 


2.19 


Kepler-18c 


0.070 ±0.107 


0.66 


Kepler-18d 


0.067 ±0.103 


0.65 



Difference image analysis takes the difference between 
average in-transit pixel images and average out-of-transit 
images. Barring pixel-level systcmatics, the pixels with 
the highest flux in the difference image will form a star 
image at the location of the transiting object, with am- 
plitude equal to the depth of the transit. Performing 
a fit of the the Ke pler pixel response function (PRF) 
(|Brvson et al.1 l2010h to both the difference and out-of- 
transit images quantifies the offset of the transit source 
from Kepler-18. Difference image analysis is vulnerable 
to various systematics due to crowding and PR F errors 
which will bias the result (jBrvson et al.1 l2010t ) . These 
types of biases will vary from quarter to quarter. We 
ameliorate these biases by computing the uncertainty- 
weighted robust average of the source locations over 
available quarters. Table Q] gives the offsets of the tran- 
sit signal source from Kepler-18 averaged over quarters 
1 through 8 for all three planet candidates. We see that 
the average offsets are within 1 sigma of Kepler-18, with 
Kepler-18b being just over 2a. From all of these lines 
of evidence, we conclude that there are no other objects 
withing 4 mags near K00137 from 0.05 arcsec out to 15 
arcsec. 

3.3. Precise Doppler Measurements of Kepler-18 

After completion of the reconnaissance spectroscopy 
and high resolution imaging, K00137 showed no evi- 
dence that might refute the planetary nature of the 
transit event in the Kepler lightcurve and the target 
was approved for high precision radial velocity obser- 
vations. We obtained 14 r elative velocities with the 
Keckl HIRES spectrometer ()Vogt et all 119941 ) between 
2009 September 1 and 2010 August 28 UT. We used 
the same spectrograph configuration as is normally used 
for precise Doppler measurements of sola r type stars 
(cf. iMarcv et alj|2008l: ICochran et al.ll2002f ). Decker B 5 
(0.87" x 3.5") was used for the first five RV observations 
taken in 2009, and decker C2 (0.87" x 14") was used for all 
of the later spectra. Exposures taken with the B5 decker 
entrance to the HIRES spectrometer suffer RV errors of 
up to ±15 ms -1 when the moon is full and the relative 
Doppler shift of the star and solar spectrum is less than 
lOkms -1 (i.e. within a line width). This 15 ms -1 er- 
ror was determined from 10 stars for which observations 
were taken with and without moonlight subtraction. The 
first 5 RV measurements of Kepler-18 were made with 
the B5 decker and hence suffer from such errors. Decker 
C2 is long enough to permit us to do accurate subtrac- 
tion of any moonlit sky spectrum. Sky subtraction was 
not possible with the spectra taken with decker B5. The 
iodine absorption cell was used as the velocity metric. 
This HIRES configuration can give a velocity precision as 
good as 1.0ms -1 , depending on the stellar spectral type, 
the stellar rotation velocity Tasini, and on the signal- 



TABLE 3 
Keck HIRES Relative Radial 
Velocity Measurements of 
K00137. 



BJD 


RV 


a 




i 

m s 


_ i 

ms 


2455076 00871 9 


4.67 


5.02 


2455076.927037 


4.24 


4.99 


2455081.024475 


5.43 


8.10 


2455082.007257 


1.09 


4.70 


2455084.983537 


-8.83 


5.61 


2455318.066020 


1.65 


5.36 


2455322.029214 


-11.97 


4.64 


2455373.003937 


10.85 


4.52 


2455403.018946 


21.47 


5.85 


2455405.909151 


-12.24 


4.73 


2455406.881186 


-0.53 


4.39 


2455413.010870 


-10.44 


5.24 


2455432.969613 


1.21 


4.45 


2455436.781954 


-8.54 


4.45 



to-noise ratio of the observation. The exposure times 
for the K00137 spectra ranged from 2300 to 2700 sec- 
onds, and the final signal-to-noise ratios ranged from 63 
to 80 per pixel, or 127 to 162 per 4.1 pixel r esolution ele- 
ment . The Doppler RV an alysis algorithm (|Butler et al.l 
119961 : Uohnson et al.1 12009| ) computes an uncertainty in 
each data point from the variance about the mean of 
the individual spectral chunks into which the spectrum 
is divided. Main sequence stars having Teff near 5400K 
have been measured for precise RVs for roughly 100 stars 
using the same HIRES instrument. Such stars typically 
show an additional noise of 2 m/s, caused by surface ve- 
locity fields and instrumental effects. We have modeled 
this additional intrinsic short-term stellar variability by 
adding that 2.0 ms -1 "jitter" in quadrature to the uncer- 
tainties of each RV data point computed by the RV code. 
The measured HIRES relative velocities are given in Ta- 
ble [1] These velocity measurements are shown as the 
black points in Figure |H The solid blue line in this fig- 
ure is the model fit from the joint MCMC solution of the 
RV data and the lightcurve, presented in Section |4j The 
rms scatter of the RV observations around this model fit 
is 4.3ms -1 , whereas the rms scatter of the raw RVs is 
9.6ms -1 . The inferred planetary masses from this solu- 
tion are given in the first row of Table [T] Also shown as 
the dashed red line in Figure |4] is a two-planet radial ve- 
locity solution that includes only Kepler-18c (K00137.01) 
and Kepler-18d (K00137.02). The black line in Figure g] 
shows the velocities from the full multi-body dynamical 
solution discussed in Section |6l The masses from this 
dynamical solution are adopted in Section [5] as our best 
determination of the planet masses. 

3.4. Spectroscopic Analysis 

We determined stellar parameters using the local ther- 
modynamic equilibrium (LTE) line analysis and spec- 
tral synthesis code MOOG 29 (jSnedenl [l973h . together 
with a grid of Kurucz 30 ATLAS9 model atmospheres. 
T he method used is virtual ly identical to that described 
in iBrugamver et al.1 (|2011l ). We analyzed a spectrum 
of Kepler-18 obtained with the Keck 1 HIRES spectro- 



available at http://www.as.utexas.edu/chris/moog.html 
http://kurucz.harvard.edu/grids.html 
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Fig. 4. — Observed precise relative radial velocity measurements 
from the Keck HIRES spectrograph. The solid black line gives 
the relative radial velocity variations of Kepler- 18 computed from 
the full multi-body dynamical model discussed in Section [6] The 
dotted (blue) line is the model fit from the combined MCMC so- 
lution to the Kepler lightcurve and these velocity data presented 
in Section [4] A Keplerian fit from the radial velocity data alone, 
without the lightcurve constraints, is indistinguishable from this 
curve. The dashed (red) line is a two-planet Keplerian fit to the 
RV data accounting for Kepler-18c (K00137.01) and Kepler-18d 
(K00137.02) only. 

graph on 27 August 2010 UT. We first measured the 
equivalent widths of a carefully-selected list of 53 neu- 
tral iron lines and 13 singly-ionized iron lines in a spec- 
trum of the Jovian satellite Ganymede, taken using the 
same instrumental setup and configuration as that used 
for Kepler-18. MOOG force-fits abundances to match 
these measured equivalent widths, using declared atomic 
line parameters. By assuming excitation equilibrium, 
we constrained the stellar temperature by eliminating 
any trends with excitation potential; assuming ionization 
equilibrium, we constrained the stellar surface gravity by 
forcing the derived iron abundance using neutral lines to 
match that of singly-ionized lines. The microturbulent 
velocity £ was constrained by eliminating any trend with 
reduced equivalent width (W\/X). Our derived stellar 
parameters for the Sun (using our Ganymede spectrum) 
are as follows: T ofr = 5785 ± 70 K, \ogg = 4.54 ± 0.09 dex, 
microturbulent velocity £ = 1.17 ± 0.06 kms" 1 , and Fe 
abundance log(e) = 7.54 ± 0.05 dex. 

We repeated the process described above for the spec- 
trum of Kepler-18. We then took the difference, on a 
line-by-line basis, of the derived iron abundance from 
each line. Our resulting iron abundance is therefore dif- 
ferential with respect to the Sun. To estimate the ro- 
tational velocity of Kepler-18, we synthesized three 5A 
wide spectral regions in the range 5640-5690A and ad- 
justed the Gaussian and rotational broadening parame- 
ters until the best fit (by eye) was found to the observed 
spectrum. The results of our analysis yield the following 
stellar parameters for Kepler-18: T c s — 5345 ± 100 K, 
\ogg = 4.31 ± 0.12, i = 1.09 ± 0.08 kms^ 1 , [Fe/H] = 



+0.19 ± 0.06 and V smi < 4kms~ 1 . The sky-projected 
rotational velocity of the star is very small. For such 
small V sini values, disentangling line-broadening effects 
due to the instrument, macroturbulence and rotation is 
difficult, at best, and requires higher signal-to-noise and 
resolution than our spectrum offers. In our MOOG anal- 
ysis, we have therefore chosen to quote an upper-limit for 
Tasini, which we estimate by assuming all broadening is 
due to stellar rotation. 

A completely independent analysis of the same spec- 
trum was done using the stellar spectral synthesis pack- 
age S ME (jValenti k Piskunovl 119961 : IValenti k Fischerl 
I2005h . This analysis of Kepler-18 gives T eff = 5383 ± 
44K, logg = 4.41 ± 0.10, [Fe/H] = +0.20 ± 0.04 and 
V smi = 0.4 ± 0.5kms _1 . The agreement between the 
MOOG and the SME analyses of this star is excellent. 
We have arbitrarily selected the MOOG parameters as 
our adopted values. 

4. LIGHT CURVE AND RADIAL VELOCITY SOLUTION 

We performed a joint solution of the photometric light 
curve and the radial velocity data in order to derive the 
orbital parameters as well as the physical characteristics 
of the planet candidates, in t he sam e manner as described 
in detail by iBatalha et all (|2011| ) for Kepler-lOb. The 
transit ligh tcurves were modeled using the analytic for- 
mulation of iMandel k Agol (|2002t ) , and the radial veloc- 
ities were fit with a Keplerian orbit. We model the mean 
density of the star, and for each of the three planet candi- 
dates we model the planetary radius, the orbital period, 
To, the impact parameter b and the radial velocity ampli- 
tude K. The eccentricity was fixed to 0.0. Model param- 
eters were determined by minimizing the \ 2 statistic us- 
ing a Levenberg-Marquardt technique. To determine the 
best-fit model parameter distributions, we used a Markov 
Chain Monte Carlo (MCMC) method th at has been o 
timized for highly correlated parameters (|Gregoryll2011l 
Stellar parameters were determined in a se parate MCMC 
analysis using the Yonsei-Yale isochrones (|Yi et al.ll200lL 
et seq.) with T c ff, [Fe/H] and fitted value of the mean- 
stellar density as constraints. Figure [5] shows the pa- 
rameter distribution functions resulting from the MCMC 
analysis. Table [1] lists the adopted system parameters. 
In most cases, these parameters came from the MCMC 
analysis. As discussed in Section [6J significant tran- 
sit timing variations were detected in Kepler-18c and 
Kepler-18d. These transit timing variations would have 
been absorbed into the MCMC periods and epochs and 
their uncertainties. Therefore we have adopted the pe- 
riods and epochs from the TTV analysis, and these are 
reported in Table [TJ The uncertainty in the epoch is 
the median absolute deviation of the transit times from 
this ephemeris and the uncertainty in the period is this 
quantity divided by the number of orbits between the 
first and last observed transits. Similarly, the dynamical 
system model presented in Section [5] derives our best and 
most reliable masses for the three planets. These masses 
from the dynamical model are reported in Table [TJ The 
masses from the MCMC model are listed in Table [TJ The 
TTV dynamical model places tight upper limits on the 
orbital eccentricity of planets " c" and " d" , as shown in 
Table [TJ These eccentricities are small enough that they 
would not significantly change the parameters in Table [TJ 
For example, the zero eccentricity solution gives a stellar 
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density of 1.01 ± 0.12 g/cm 3 . Inclusion of the eccentricity 
from the TTV solution would change this by 0.03 g/cm 3 , 
or 0.25cr. For each parameter adopted from the MCMC 
model, we adopt the median value of the distribution, 
and the error bar represents the 68.3% confidence level 
- roughly equivalent to a lcr confidence level. 

It is interesting and informative to compare the initial 
values for the transiting planet p arameters pres e nted in 
the tabulation of candidates by iBorucki et al.1 (|2011b[ ) 
with those we have finally settled on in Table Q] The 
parame ters of the parent star given by IBorucki et al.l 
(|2011b[ ) K00137 are simply taken from the KIC. The 
overall accuracy of the KIC photometric st ellar parame- 
ters is actually quite good, as was shown bv lBrown et al.l 
(|2011f) . The othe r important p l anetary system param- 
eters reported by IBorucki et al.l ()2011bl ) are the orbital 
periods and transit depths of each of the three KOIs in 
the system. Our adopted o rbital p eriods agree with the 
preliminary IBorucki et al.l ()2011bf ) values to within 1 a 
for K00137.01 and K00 137.03, and w i thin 4 sigma for 
K00137.02. We note that IBorucki et al.l (|2011bl) assumed 
no transit time variations, and fit the best mean period 
to the data. The transit depths agree to within 2 a for 
K00137.02 and K00137.03, but differ by 7.5 sigma for 
K00137.02. This slight disagreement fo r K00137.02 prob- 
ably results from IBorucki et all (|2011bf ) fitting an impact 
parameter of 0.03, whereas our best value is 0.593±0.050. 

The radial velocity variations of the parent star Kepler- 
18 resulting from this three planet MCMC model is 
shown as the solid blue line in Figure |U This solution 
is very similar to the three-body radial velocity solution 
one would get from the radial velocity data alone, hold- 
ing the period and T of each planet fixed at the values 
from the photometric lightcurve solution and assuming 
zero eccentricity. A two-body solution to the radial ve- 
locity data considering only K00137.01 and K00137.02 
is shown as the dashed red line in Figure @] The full 
three-planet solution gives a significantly better fit to 
the observed radial velocities than does the two-body 
solution. The rms of the full MCMC three-body fit is 
4.3 ms -1 , while the rms of the two-body Keplerian solu- 
tion is 5.5ms -1 . We computed an /-test based on the 
residuals of a 2 planet and 3-planet RV fit. The periods 
where fixed based on transit data. We find / = 4.97, 
which means there is a 0.68% probability that variance 
of the residuals are similar. Thus this reduction in the 
rms supports the interpretation of K00137.03 as a third 
planetary companion to K00137. However, the radial ve- 
locity observations are not sufficiently dense to consider 
this to be a confirmation of this KOI as a planet. 

5. Warm-Spitzer OBSERVATIONS OF K00137.01 AND 
K00137.02 

Observation of the transits of the planets around 
Kepler-18 at multiple, widely-spaced wavelengths is a 
valuable tool to confirm the planetary nature of the 
events. The depth of a planetary transit should be 
nearly independent of wavelength, aside from minor ef- 
fects due to the possible finite brightness of the planet 
as a function of wavelength. In order to test for wave- 
length independence of the transit depths, K00137.01 
and K00137.02 were o bserved during two t ransits with 
Warm -Spitzer /IRAC (fWerner et al.l 12004 iFazio et al.l 
120041) at 4.5 /^m (program ID 60028). The observa- 



TABLE 4 

Adopted System Parameters 



Parameter 



Value 



(M ) a 

r* (R@r 


0.972 ± 0.042 


1.108 ± 0.051 


log!,* (L ) a 
Age (Gyr) a 


-0.031 ± 0.035 


10.0 ± 2.3 


logg* b 


4.31 ±0.12 


P* (g/cm 3 ) c 


1.01 ± 0.12 



Kepler-18b 

T C 

P (days) c 

Transit depth (ppm) c 
b (Impact Parameter) 
Rp/ R+ c 

M p (M ffi ) d 

Rp (R®) c 

i (dcg)° 

Pp (g/cm 3 ) 3 ' 4 
a (AU) C 
K (m/s) c 
Tw (h) e 
T12 (h) c 

Kepler-18c : 

T d 

P (days) d 

Transit depth (ppm) c 
b (Impact Parameter) c 
Rp/R* c 
M p (M ffi ) d 

Rp (R(bT 

% (deg)° 
a/R* e 

Pp (g/cm 3 ) 3 ' 4 
a (AU) C 
K (m/s) c 
Tw (h) e 
T12 (h) c 



K00137.03 
2454966.5068 ± 0.0021 
3.504725 ± 0.000028 
254.0 ±7.8 
0.771 ± 0.025 
0.01656 ±0.00032 
6.9 ±3.4 
2.00 ±0.10 
84.92 ± 0.26 
8.58 ±0.37 
4.9 ±2.4 
0.0447 ± 0.0006 
5.2 ±2.4 
2.076 ± 0.036 
0.0818 ±0.0082 
: K00137.01 
2455167.0883 ± 0.0023 
7.64159 ± 0.00003 
2286.6 ±8.6 
0.593 ±0.050 
0.04549 ± 0.00055 
17.3 ± 1.9 
5.49 ± 0.26 
87.68 ± 0.22 
14.43 ± 0.61 
0.59 ±0.07 
0.0752 ±0.0011 
5.1 ± 1.9 
3.488 ±0.020 
0.229 ±0.022 



Kepler-18d 
T d 

P (days) d 

Transit depth (ppm) c 
b (Impact Parameter) c 
R P /R* C 
M P (M e ) d 

Rp (R®r 

i (dcg)° 
a/R* e 

Pp (g/cm 3 ) 3 ' 4 
a (AU) C 
K (m/s) c 
T14 (h) e 
Tl2 (h) e 



= K00137.02 
2455169.1776 ±0.0013 
14.85888 ± 0.00004 
3265. ± 12. 
0.767 ±0.024 
0.05782 ± 0.00069 
16.4 ± 1.4 
6.98 ±0.33 
88.07 ±0.10 
22.48 ± 0.96 
0.27 ±0.03 
0.1172 ±0.0017 
7.3 ±2.1 
3.679 ± 0.036 
0.459 ± 0.045 



a based on isochrone fits using p t , from MCMC 

model and T e ff and [Fe/H] from spectroscopy. 

b from MOOG analysis 

c from MCMC model 

d from TTV model 

e derived from other parameters 

tions occurred on UT 2010 July 19 and UT 2010 August 
13. Both visits lasted approximately 9h. The data were 
gathered in full-frame mode (256 x 256 pixels) with an 
exposure time of 12 s per image which yielded 2418 and 
2575 images per respective visit. 

The method we used to produce photometric t i me se - 
ries from the images is described by Desert et al. ( 2009) . 
It consists of finding the centroid position of the stel- 
lar point spread function (PSF) and performing aperture 
photometry using a circular aperture 

The images used are the Basic Calibrated Data (BCD) 
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Fig. 5. — Markov Chain Monte Carlo model probability distribution functions for parameters of the Kepler-18 system. The left column 
gives the two stellar parameters of the stellar density and the relative radial velocity zero point. For the other five columns, the top row 
gives the planet parameters of transit epoch, orbital period, square of the impact parameter, ratio of the planetary to stellar radius, and 
the orbital velocity amplitude K for Kepler-18b. The second and third rows give the same parameters for Kepler-18c and Kepler-18d 
respectively. 

dures, GCNTRD, from the IDL Astronomy Library 31 . We 
use the APER routine to perform aperture photometry 
with a circular aperture of variable radius, using radii 
of 1.5 to 8 pixels, in 0.5 pixel steps. The propagated 
uncertainties are derived as a function of the aperture 
radius; we adopt the one which provides the smallest er- 
rors. We find that the transit depths and errors vary only 
weakly with the aperture radius for all the light-curves 
analyzed in this project. The optimal apertures is found 
to be at 3.0 pixels. We estimate the background by fit- 
ting a Gaussian to the central region of the histogram 
of counts from the full array. The center of the Gaus- 
sian fit is adopted as the residual background intensity. 
As already seen in previous Warm- Spi tzer observations 
(|Deming et al.ll201ltlBeerer et alj|2011l) . we find that the 
background varies by 20% between three distinct levels 
from image to image, and displays a ramp-like behav- 
ior as function of time. The contribution of the back- 
ground to the total flux from the star is low for both 
observations, from 0.15% to 0.9% depending on the im- 
ages. Therefore, photometric errors are not dominated 

31 http : // idlastro . gsf c. nasa . gov/homepage . html 



delivered by the Spitzer archive. These files are corrected 
for dark current, flat-fielding, detector non-linearity and 
converted into flux units. We convert the pixel intensities 
to electrons using the information given in the detector 
gain and exposure time provided in the FITS headers. 
This facilitates the evaluation of the photometric errors. 
We extract the UTC-based Julian date for each image 
from the FITS header (keyword DATE.OBS) and correct 
to mid-exposure. We convert to UTC-based BJD follow- 
ing the procedure developed by lEastman et all {2010). 
We use the JPL Horizons ephemeris to estimate the po- 
sition of Spitzer Space Telescope during the observations. 
We correct for transient pixels in each individual image 
using a 20-point sliding median filter of the pixel inten- 
sity versus time. For this step, we compare each pixel's 
intensity to the median of the 10 preceding and 10 fol- 
lowing exposures at the same pixel position and we re- 
place outliers greater than 4a with their median value. 
The fraction of pixels we correct is lower than 0.16% for 
both transits. The centroid position of the stellar PSF 
is determined using DAOPHOT-type Photometry Proce- 
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by fluctuations in the background. We used a sliding me- 
dian filter to select and trim outliers in flux and position 
greater than 5 cr, which corresponds to 0.9% and 2.0% of 
the data, for the first and second visit respectively. We 
also discarded the first half-hour of observations, which 
are affected by a significant telescope jitter before stabi- 
lization. The final number of photometric measurements 
used is 2246 and 2369. The raw time series are presented 
in the top panels of Figure [51 We find that the typical 
signal-to-noise ratio (S/N) 140 per image which corre- 
sponds to 80% of the theoretical signal-to-noise. There- 
fore, the noise is dominated by Poisson photon noise. 

5.1. Analysis of the Warm-Spitzer light curves 

We used a transit light curve model multiplied by in- 
strumental decorrelation functions to measure the tran- 
sit parameters and th eir uncertain ti es from the Spitzer 
data as described in IDesert et al.l (|2011bl) . We com- 
pute the transit light c urves with the IDL tra nsit rou- 
tine OCCULTSMALL from iMandel fc Agoll (|2002h . In the 
present case, this function depends on one parameter: 
the planet-to-star radius ratio R p / R+. The orbital semi- 
major axis to stellar radius ratio (system scale) a/R+, the 
impact parameter b, and the time of mid transit To are 
fixed to the values derived from the Kepler lightcurves. 
The limb-darkening coefficients are set to zero since these 
Spitzer lightcurves do not have enough photometric pre- 
cision to detect the curvature in the transit light curve 
that would be produced by limb darkening. 

The Spitzer /IRAC photometry is known to be sys- 
tematica lly affected by the so-ca l led pixel-phase effec t 
(see e.g.. lCharbonneau et al.ll2005t iKnutson et al.1 12008) . 
This effect is seen as oscillations in the measured fluxes 
with a period of approximately 70 min (period of the tele- 
scope pointing jitter) and an amplitude of approximately 
2% peak-to-peak. We decorrelated our signal in each 
channel using a linear function of time for the baseline 
(two parameters) and a quadratic function of the PSF po- 
sition (four parameters) to correct the data for each chan- 
nel. We checked that adding paramet ers to the correction 
function of the PSF position (as in De sert et al.ll2009| ) 
does not improve the fit significantly. We performed 
a simultaneous L evenberg-Marquardt least-squares fit 
(jMarkwardtl 1200 9) to the data to determine the tran- 
sit and instrumental model parameters (7 in total). The 
errors on each photometric point were assumed to be 
identical, and were set to the rms of the residuals of 
the initial best-fit obtained. To obta in an estimate o f 
the correlated and systematic errors (|Pont et all l2006h 
in our measurements, we use the residual permutation 
bootstrap, o r "Pra yer Bead", method as described in 
IDesert et al.l ((2009) . In this method, the residuals of 
the initial fit are shifted systematically and sequentially 
by one frame, and then added to the transit light curve 
model before fitting again. We allow asymmetric error 
bars spanning 34% of the points above and below the me- 
dian of the distributions to derive the 1 a uncertaintie s 
for each parameters as described in lDesert et al.1 ()2011al ). 

We measured the transit depths at 4.5 ftm of 152ll245 
ppm for Kepler-18c and 40831^85 for Kepler-18d. The 
values we measure for the transit depths in the Spitzer 
bandpass are in agreement at the 2a level compared to 
the Kepler bandpass. This indicates that the transit 



depths of Kepler- 18c and Kepler- 18d are only weakly de- 
pendent on wavelength, if at all. This is in agreement 
with expectations for a dark planetary object, and indi- 
cates that there is no significant contamination from any 
nearby unresolved star of significantly different color. As 
discussed in Section[7l this is important evidence that the 
transit signals arise from planets and not from eclipsing 
stars blended with additional light. 

These Spitzer observations provide a useful constraint 
on the kinds of false positives (blends) that may be 
mimicking the transit signal, such as eclipsing binaries 
blended with the target. If Kepler-18 were blended with 
an unresolved eclipsing binary of later spectral type that 
manages to reproduce the transit depth in the Kepler 
passband, the predicted depth at 4.5/mi would be ex- 
pected to be larger because of the higher flux of the con- 
taminant at longer wavelengths compared to Kepler-18. 
Since the transit depth we measure in the near infrared 
is about the same as in the optical, this argues against 
blends composed of stars of much later spectral type. 
Based on model isochrones, the properties of the target 
star, and the transit depths measured with Spitzer at the 
3(7 level, we determine a lower limit to the blend masses 
of 0.86 and 0.79M©, for Kepler-18c and Kepler-18d re- 
spectively. 

6. TRANSIT TIMING VARIATIONS ANALYSIS 
6.1. Transit Times and Errors 

Kepler-18 was observed at the 29.4244 minute Long 
Cadence (LC) rate for the first three observing periods 
Q0-Q2. After the transits of K00137.01 and K00137.02 
were detected, Kepler-18 was observed at the 58.85 sec- 
ond Short Cadence in order to facilitate the measurement 
of possible transit timing variations. These SC data were 
used for Q3 through Q7. Transit times and their errors 
for each member of the Kepler-18 system were deter- 
mined through an iterative procedure; a single step of 
this procedure is described as follows. 

First, the detrended photometric data within four tran- 
sit durations at each epoch were shifted by the current 
best-fit mid-transit times, and the lightcurves were folded 
on these transit times to form a transit template. This 
template was then fit with a transit light curve model 
(|Mandel fc Agolll2002l) . Next, at each individual epoch, 
this light curve template was shifted in time and com- 
pared to the data by computing the standard \ 2 statistic. 
This statistic was computed for a dense, uniform sample 
of mid-transit times centered on the current best-fit time 
and spanning approximately five current best-fit timing 
errors. The time to, corresponding to the minimum % 2 , 
was recorded as the new best-fit mid-transit time. The 
X 2 data were then fit with a quadratic function of time, 
X 2 (t) — C(t — i ) 2 + Xo- By choosing this functional 
form, we are assuming that the posterior likelihood is 
well-described by a Gaussian function of the mid-transit 
time. The fidelity of the quadratic approximation was 
verified visually at each epoch. The timing error, cr t , was 
found by solving for the time, t = to + a t , at which the 
quadratic model indicated a A% 2 = 1. This corresponds 
to a t = C- 1 / 2 . 

This iterative procedure of template generation fol- 
lowed by timing estimate converged to the final values 
in generally two to three steps. The measured transit 
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Fig. 6. — Spitzer transit light-curves of K00137.01 (top) and K00137.02 (bottom) observed in the IRAC band-pass at 4.5 )im. Top 
panels : raw (unbinned) transit light-curves. The red solid lines correspond to the best fit models which include the time and position 
instrumental decorrelations as well as the model for the planetary transit (see details in Sect. [SJ. Lower panels : transit light-curve corrected, 
normalized and binned by 36 minutes. The best-fit Spitzer transit curves are plotted in red and the transit shapes expected from the 
Kepler observations are overplotted in dashed green lines. 



times for each observed transit of Kepler-18c are given 
in Table[TJ and the measured transit times for Kepler-18d 
are given in Table [T] 

A drift in the orbital inclination was measured for each 
planet by augmenting the nominal transit model with an 
epoch-dependent linear inclination (viz., i(E) — i(0) + 
(Ai/AE) x E), and then fitting for the linear coefficient. 
In this fit, transit times were fixed to their best- fit values 
while the remaining transit parameters were allowed to 
vary. The best-fitting solution was found by minimizing 
the standard % 2 metric. The uncertainty was estimated 
by fitting a multivariate Gaussian to a sampling of the 
posterior parameter distribution. The drifts (Ai/AE) 
were found to be [10 ± 9, 4 ± 22, -3 ± 10] x 10~ 4 degrees 
per epoch for planets b, c, and d respectively. Thus we 
detect no secular inclination drift. 

6.2. Transit Time Variation Analysis 

To model the transit times and radial velocities, we 
applied the numerical routines that were previously used 
for Kepler-9 and -11 pa"brvcky||2010l:lHolman et al'.ll2?J10r, 
iLissauer et "all 1201 lal) . That is, a Levenberg-Marquardt 
^-minimization routine drove 3-planet dynamical in- 
tegrations, which calculated the transit times resulting 
from g iven orbital parameters. As in ILissauer et al.l 
(|2011al ). we chose the parameters (m p , P, T , ecosw, 
esinw) for each planet: mass, orbital period, transit 
phase, and eccentricity vector components. These pa- 
rameters are osculating Jacobi an orbital elements a t the 
epoch 2455168.0 [BJD]. As in iHolman et al.l (|2010D . we 
also used the integrations to calculate radial velocities 
of the star at each of the observed times. We assumed 
M* = 0.95M Q , assumed the orbits are edge-on and copla- 
nar, and neglected light travel time effects in these nu- 
merical calculations. 

The main data constraining the orbital model are the 
transit times of Tables Q] (75 data points) and [T] (37 
data points). In some of the calculations reported be- 
low, we allow K00137.03 to interact dynamically, but 
we include only its first and last observed transits, at 
t = 2454955.99237±0.00823and 2455530. 77178±0.00174 



(2 data points) , as a way of keeping its period and phase 
fixed at observed values. We have measured all of its 
transit times, but we do not report or analyze them here: 
their individual signal to noise is low, they are not ap- 
parently constant, and we have not found a consistent 
dynamical solution to date. Perhaps future work will 
show a fourth (non-transiting) planet is required to fit 
the transit times of Kepler- 18b. In the meantime we 
proceed with fitting the three transiting planets with a 
focus on the transit timing constraints for Kepler-18c and 
Kepler- 18d. 

By fitting only the transit times, allowing P and To 
for each of the three planets to vary (6 free parame- 
ters), and setting dynamical interactions to zero, we find 
X 2 = 750.6 for 114 transit time measurements. This is 
clearly an unacceptable fit, and the obvious timing pat- 
terns call for a dynamical model. The observed ("O") 
residuals to this calculated ("C") ephemeris are called 
the O-C values, and are plotted in Figure [7] along with 
the preferred dynamical model described below. As in 
the case of Kepler-llb/c, we clearly see that the source of 
the variations of transit times for Kepler-18c and Kepler- 
18d is their near 2:1 resonance. The expected variation 
occurs on a timescale: 

Pttv = l/(2/Pd - l/Pc) = 268 days, (1) 

the time it takes the line of conjunctions to sweep around 
inertial space, through both the line of sight and the apsi- 
dal lines of these planets (on the approximati on that pre- 
cessio n can be ignored on this timescale); see lAgol et al.l 
( 2005). In figure [U we plot the periodogram of the O-C 
values for Kepler-18c and Kepler-18d. The peaks oc- 
cur at 1/(260.4 ± 3.3 days) and 1/(265.1 ± 4.7 days) for 
Kepler-18c and Kepler-18d respectively, very close to the 
simple expectation given above, which uniquely identifies 
the dynamical mechanism for transit timing variations. 

Before moving on to a full solution, we tried fitting the 
radial velocities (14 data points, Tabled]) with the above- 
determined periods and phases, with circular orbits, al- 
lowing only the planetary masses to vary. The solution, 
listed in the second row of Table [I] was [m?,, m c , m^] = 



Kepler-18b,c,d 



11 



TABLE 5 

Transit Times for Kepler-18c = K00137.01 



cycle a 


BJD-2455000.0 


O-C 


<7 


-27.0 


-39 


.23132 


-i-o oo^4 ,; i 


001 41 


-26.0 


-31 


.59238 


_i_n 00077 


001 28 

KJ • VI KJ liu 


-25.0 


-23 


.94944 


4-0 0091 9 


001 14 

U.UU-LO'L 


-24.0 


-16 


.30829 


4-0 001 fi8 


001 9^ 

U.UU-LZU 


-23.0 


-8 


.66703 


4-0 001 14 

I KJ .UU1 'Jt: 


00909 

KJ .Uv^Ui 


-21.0 


6 


.61139 


— 001/19 

— U.UUO'±^ 


001 fi1 

U . UU -LU -L 


-20.0 


14 


.25400 


_n 00941 


001 ^fi 

U. UU-LUU 


-19.0 


21 


.89543 


— 009^7 


001 4Q 

U. UU-L^if 


-18.0 


29 


.53633 


— 00197 


001 1 1 ^ 

U . UU -LOU 


-17.0 


37 


.17764 


_n 001^ 

KJ . VJKJiJ^J ij 


001 11 

U. UU1U -L 


-16.0 


44 


.82280 


4-0 00009 


001 A(\ 

U. UU-L^iU 


-15.0 


52 


.46001 


— 0041fS 


001 07 

KJ . VI KJ _LU 1 


-14.0 


60 


.10428 


— 001 fiQ 


001 1 ^ 

U . UU -L J. O 


-13.0 


67 


.74533 


— 00991 


001 08 

U.UU1UO 


-12.0 


75 


.39019 


4-0 001 01 

1 U .UU1UU 


001 4^ 

KJ . VI KJ A-'-t'J 


-11.0 


83 


.02837 


— 00918 


001 44 

U. UU14:4 


-10.0 


90 


.67043 


— 001 Q1 

U . UU -L u -L 


001 40 

U. 


-9.0 


98 


.31334 


— OOO^Q 


001 09 

U. UU1UZ 


-8.0 


105 


.95509 


— 00044 


001 O^i 

U . UU -LUU 


-7.0 


113 


.59742 


4-0 00010 

1 KJ . WUu JU 


001 48 

u. WlW 


-6.0 


121 


.23989 


4-0 001 1 7 

TU.UUl J. t 


000Q8 
u.uuuyo 


-5.0 


128 


.88188 


4-0 001 ^7 

1 KJ . UU _LU 1 


001 08 

U . UU -L UU 


-4.0 


136 


.52252 


4-0 000R9 


001 09 

U . UU -LU^ 


-3.0 


144 


.16670 


4-0 00190 

|^ KJ . KJKJ'J i-J\J 


00088 

U. uUUUU 


-2.0 


151 


.80782 


4-0 00273 

\^ KJ . UUii 1 t > 


001 37 

U . UU 1J 1 


-1.0 


159 


.44986 


+0 0031 8 


001 1 3 

U . UU -L J-O 


0.0 


167 


.09286 


+0 001 59 

1 KJ . KJKJ'-t'JZf 


00096 

U. UUU i/U 


1.0 


174 


.73457 


4-0 00470 

TU.UUl 1 U 


000Q7 
u.uuuy I 


2.0 


182 


.37715 


4-0 OO'SfSX 


001 11 

U . UU lul 


3.0 


190 


.01671 


+0 00365 

1 KJ . UUUU 


001 06 

U. UU1UU 


4.0 


197 


.65742 


4-0 00277 


001 09 

U. UU-LUC/ 


5.0 


205 


.30025 


4-0 00400 


0.00117 


6 


212 


.94064 


4-0 00981 


001 01 

U . UU -LU -L 


7.0 


220 


.57949 


4-0 OOOOfi 

I KJ . UUUUU 


0.00111 


8.0 


228 


.22281 


4-0 001 78 

TU.UUl I O 


001 01 

U.UU1UL 


9.0 


235 


.86365 


4-0 001 01 

1 KJ ,UU1UU 


001 OfS 

U . UU 1UU 


10.0 


243 


.50537 


4-0 001 1 

TU.UUl -LU 


001 07 

U . UU -LU t 


11.0 


251 


.14597 


4-0 0001 7 

TU.UUUl t 


000Q7 
u.uuuy I 


12.0 


258 


.78515 


_0 00994 


000Q9 

U . VI K J VI CJ 


13.0 


266 


.42783 


— 001 1 

— U.UU11U 


001 1 ^ 

U . UU -L J- O 


14.0 


274 


.06799 


— 009^9 

KJ . KJKJA'JCf 


001 1 8 

U. uUl lu 


15.0 


281 


.70882 


— 0011^ 

— U. UUOOtJ 


001 1 

U.UU11U 


16.0 


289 


.35100 


— 00977 


0007Q 
u . uuu i y 


17.0 


296 


.99249 


— 0098X 


00088 

U. uUUUU 


18.0 


304 


.63304 


— 001Q9 

— u.uuoy^ 


00081 
u. uuuoo 


19.0 


312 


.27505 


— 001^0 

KJ . UUOUU 


000Q1 

U . VIKJV) CJ ' J 


20.0 


319 


.91502 


— 00^*1 9 

— U.UUJli 


000Q0 
u . uuuyu 


21.0 


327 


.56008 


— 001 F,^ 

— U.UU1UO 


0008Q 
u.uuuoy 


22.0 


335 


.19905 


— 00497 


000Q1 

U . VIKJV) CJ 'J 


23.0 


342 


.84138 


— 001^ 

— U. UUOUO 


00084 
u.uuuoi 


24.0 


350 


.48436 


— 0091 5 

KJ . UUi J- <J 


001 01 

U. UUi KJ'J 


25.0 


358 


.12585 


— 0099R 


000Q7 
u.uuuy I 


26.0 


365 


.76604 


— 001R7 

— U.UUOU t 


OOOSfi 
u. uuuou 


27.0 


373 


.41078 


_0 000^9 

KJ . UUUUi 


OOOQ'S 

U. UUU Ju 


28.0 


381 


.05287 


— 00001 

— U.UUUU1 


001 01 

U.UU1UL 


29.0 


388 


.69546 


4-0 000Q8 

\^ KJ . KJKJVJTyKJ 


00089 

U. KJKJVK3 U 


30.0 


396 


.33836 


4-0 00998 

TU.UUZ^O 


000Q0 
u.uuuyu 


31.0 


403 


.97897 


4-0 001 10 

TU . UU -LOU 


OOOQfi 
u.uuuyu 


32.0 


411 


.62239 


4-0 0011 1 


001 04 

U . UU lUi 


33.0 


419 


.26374 


4-0 00988 
TU . uuzoo 


000Q0 
u.uuuyu 


34.0 


426 


.90502 


4-0 009^6 

1 KJ .WUiiJU 


001 00 

U. UUi UU 


35.0 


434 


.54873 


4-0 004R8 

TU.UU1UO 


001 1 9 

U. UUI LZi 


36.0 


442 


.18944 


4-0 00181 

TU.UUOO -L 


000Q1 
u.uuuy -L 


37.0 


449 


.82907 


4-0 001 84 

I KJ .UUlOi 


001 OfS 

U. UUi UU 


38.0 


457 


.46935 


+0.00052 


0.00097 


39.0 


465 


.11288 


+0.00246 


0.00092 


40.0 


472 


.75469 


+0.00268 


0.00103 


41.0 


480 


.39522 


+0.00161 


0.00117 


42.0 


488 


.03539 


+0.00020 


0.00101 


43.0 


495 


.67681 


+0.00002 


0.00100 


44.0 


503 


.31876 


+0.00037 


0.00092 


45.0 


510 


.95872 


-0.00125 


0.00091 


46.0 


518 


.60144 


-0.00014 


0.00088 


47.0 


526 


.24067 


-0.00249 


0.00114 


48.0 


533 


.88235 


-0.00241 


0.00112 



TABLE 6 

Transit Times for Kepler-18d = K00137.02 



cycle a BJD-2455000.0 O-C 



-14 


.0 


-38 


O A f\ CI f\ 

.84989 


— 0. 


.00318 


0. 


.00101 


-13 


.0 


-23 


.98927 


— 0. 


.00144 


0. 


.00090 


-12. 


.0 


-9 


.12888 


+0. 


.00008 


0. 


.00081 


-11 


.0 


5 


.73040 


+0 


.00048 


0. 


00086 


-10 


.0 


20 


.59242 


+0 


r\ r\ o c o 

.00362 


0. 


00115 


-9 


.0 


35 


A A OC\ A 

.44894 


+0 


.00126 


0. 


00074 


-8 


.0 


50 


.30774 


+0 


.00119 


0. 


00091 


-7 


.0 


65 


.16644 


+0 


.00101 


0. 


00078 


-6 


.0 


80 


.02557 


+0 


.00126 


0. 


00121 


-5 


.0 


94 


.88445 


+0. 


.00127 


0. 


U0U68 


- 1 


.0 


109 


.74201 


— 0. 


.00005 


0. 


.00065 


-3 


.0 


124 


.59986 


— 0. 


.00107 


0. 


.00080 


-2 


.0 


139 


.45743 


— 0. 


f\ c\ o o o 

.00238 


0. 


.00065 


-1 


.0 


154 


.31562 


— 0. 


f\f\ OA *~7 

.00307 


0. 


.00067 





.0 


169 


.17495 


— 0. 


.00261 


0. 


U0U7U 


2 


.0 


198 


.89361 


— 0. 


.00170 


0. 


.00081 


3 


.0 


213 


.75079 


— 0. 


.00341 


0. 


.00085 


1 


.0 


228 


.61267 


— 0. 


.00040 


0. 


.00079 


5. 


.0 


243 


.47039 


— 0. 


.00155 


0. 


.00081 


6 


.0 


258 


o o n o o 

.33093 


+0 


.00011 


0. 


U00d8 


7 


.0 


273 


.19074 


+0 


.00104 


0. 


00071 


8. 


.0 


288 


.05139 


+0 


.00282 


0. 


00069 


9 


.0 


302 


.91049 


+0 


.00304 


0. 


00061 


10. 





317 


.77010 


+ 


.OUo/ I 


0. 


.U0U65 


11. 





332 


.62850 


+0. 


.00330 


0. 


.00061 


12. 





347 


.48637 


+0. 


.00229 


0. 


.00069 


13. 





362 


.34423 


+0 


.00127 


0. 


.00072 


14. 





377 


.20286 


+0 


.00103 


0. 


.00070 


15. 





392 


.06096 


+0. 


.00025 


0. 


.00063 


16. 





406 


.91891 


-0. 


.00068 


0. 


00061 


17. 





421 


.77755 


-0. 


.00091 


0. 


00061 


18. 





436 


.63489 


-0 


.00245 


0. 


00063 


19. 





451 


.49447 


-0 


.00175 


0. 


00067 


20. 





466 


.35351 


-0. 


.00158 


0. 


00063 


21. 





481 


.21293 


-0 


.00104 


0. 


00066 


22. 





496 


.07232 


-0. 


.00052 


0. 


00067 


23. 





510 


.93180 


+0 


.00008 


0. 


.00061 



a P = 14.85888 days, T = 2455169.17756 

[12 ±5, 15 ± 5, 28 ± 7]M ffi ; the X 2 = 9.7 for 10 de- 
grees of freedom (14 radial velocity data points, minus 
3 K-amplitudcs, minus a constant radial- velocity offset). 
These masses are well within the la error bars of the 
MCMC solution to the lightcurve and RVs presented in 
Section [4] Therefore this model is sufficient to explain 
the radial velocities, and each planet is detected, but 
only marginally so for planet b. In particular, if we hold 
the mass of b at zero, the masses of the others become 
[m C) m d ] = [18 ± 5, 24 ± 7] M e and the \ 2 = 15.6 for 10 
degrees of freedom. These values are within about 0.5a 
of the masses determined from the joint MCMC and RV 
solution presented earlier in Section |4j 

We may also attempt to constrain the masses and or- 
bital elements using only the transit times. Naturally, 
this requires full dynamical integrations, in which the 
planets cannot remain on circular orbits. However, for 
planet b we assumed a circular orbit at the dynami- 
cal epoch, since we are not attempting to fit its transit 
time variations. All of the other orbital parameters and 
masses were free to vary. The resulting x 2 is 88.5 for 101 
degrees of freedom (114 transit times, minus 5 parame- 
ters for planets c and d, minus 3 for K00137.03), which 
is quite acceptable. To be compared with the radial- 
velocity solution, the solved- for masses were [rrib, m Cl rrid] 
= [18 ± 9, 17.3 ± 1.7, 15.8 ± 1.3]M ffi , shown in the third 
row of Table Q] These solutions had extremely low ec- 
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centricities (e < 0.003) for Kepler-18c and Kepler-18d. 
As above, the innermost planet is only very marginally 
detected. In contrast to the radial- velocity solution, how- 
ever, the masses of the interacting planets are very pre- 
cisely pinned down by the large variations seen in Fig- 
ure GO 

Finally, we generated a joint solution to the transit 
times and radial velocities. Graphically, this solution is 
given by Figure [7J Ax 2 — 103.4 for 114 degrees of free- 
dom is achieved, an excellent fit to the data. The orbital 
parameters and their formal errors (the output of the 
Levenberg-Marquardt algorithm) are given in Table [IJ 
and the planetary masses are given in the fourth row of 
Tabic CO 

In the previous subsection we found no drift in the 
inclinations of the three planets. We can use the 3cr up- 
per limits on \Ai/AE\ of [37, 70, 33]xl0~ 4 degrees per 
epoch for Kepler- 18 b, c, and d to place l imits on their 
mutual inclination (|Miralda-Escudel f2002f ) . To do this, 
we measure the value of |Ai/A_E| seen in numerical sim- 
ulation, which depends nearly linearly on the difference 
in no dal angle on the sky of two planets ()Ballard et al.l 
2010). Taking the masses and orbits from the best-fit 
TTV/RV solution above, we simulated Kepler-18c and d 
with 1 degree (and 10 degrees) of mutual inclination, we 
find \ Ai/AE\ c = 1.8 x 10~ 4 degrees per epoch (16 x 10~ 4 
degrees per epoch) and \Ai/AE\d — 2.9 x 10~ 4 degrees 
per epoch (26 x 10 -4 degrees per epoch). The interaction 
between Kepler-18b and c is also potentially observable; 
simulating their orbits with 1 (10) degree(s) of mutual in- 
clination, we find | Ai/Ai?|(,=1.3(ll) x 10~ 4 degrees per 
epoch and |Ai/A_E| c =0.82(7.4) x 10~ 4 degrees per epoch. 
The interaction between Kepler-18b and d is consider- 
ably weaker, with a drift of ~ 2 x 10 -4 degrees per epoch 
even for 10 degree mutual inclination, therefore we ignore 
it when setting mutual inclination limits. To find the 3cr 
upper limit to the nodal difference of c and d, we compare 
the 10-degree calculated value of \Ai/AE\d to its obser- 
vational 3er upper limit, so we find |f2 c — Q^l < 13°. The 
inclination difference in the complementary direction is 
only id — i c = 0.40 ± 0.24 degrees, so the limit on the 
true mutual inclination is i c d < 13 degrees. The nodal 
difference of b and c is more poorly constrained, as the 
calculated value of \Ai/AE\ for both planets in the 10- 
degree mutual inclination simulation is small compared 
to its observational 3er upper limit. Thus a moderate 
(~ 20°) mutual inclination is permissible, which would 
be large enough to affect the timing fits. In that case, 
the model should self-consistently fit the radial velocities, 
transit times, and transit durations. However, we defer 
dynamical interpretation of the mutual inclination of the 
inner planet to the others until more data are gathered, 
which should tighten the limit. 

7. BLENDER ANALYSIS OF K00137.03 

The lack of a clear dynamical confirmation of the na- 
ture of K00137.03 requires us to examine the wide variety 
of astrophysical false positives (blends) that might mimic 
the photometric transit, and to assess their a priori like- 
lihood compared to that of a true planet. For this we 
apply the BLENDER technique described by iTorres et al.l 
(2004 . 1201 If), with further developments as reported by 
iFressin et atl (|201lH . BLENDER uses the detailed shape of 
the transit light curve to weed out scenarios that lead to 



the wrong shape for a transit. The kinds of false positives 
we are concerned with here include background or fore- 
ground eclipsing binaries blended with the target, as well 
as physically associated eclipsing binaries, which gener- 
ally cannot be resolved in high-angular resolution imag- 
ing. In each case the pair of eclipsing objects can also 
be a star transited by a larger planet. Briefly, BLENDER 
simulates a very large number of light curves resulting 
from these blend scenarios with a range of stellar (or 
planetary) parameters, and compares them to the Ke- 
pler photometry in a \ 2 sense. Blends providing poor 
fits are considered to be ruled out, enabling us to place 
constraints on the kinds of objects composing the eclips- 
ing pair that yield viable blends, including their size or 
mass, as well as other properties of the blend such as the 
overall brightness and color, and even the eccentricities 
(e) of the orbits. We refer the reader to the above ref- 
erences for details. Following the nomenclature in those 
sources, the objects in the eclipsing pair are designated 
the 'secondary' and 'tertiary', and the target itself is 
the 'primary'. Stellar properties are drawn from model 
isochrones. 

Simulations with BLENDER indicate that background 
eclipsing binaries with two stellar components can only 
produce viable false positives if they are restricted to 
a narrow range of masses for the secondaries (approxi- 
mately 0.8 < M2/MQ < 1.3), as well as a limited inter- 
val in brightness (Kp magnitude) relative to the target 
(4 < AKp < 7). This is illustrated in Figure M in which 
we show the \ 2 landscape from BLENDER for all blend fits 
of this kind. Regions outside the 3-er contour correspond 
to scenarios with transit shapes that do not provide ac- 
ceptable fits to the Kepler photometry, i.e., fits that are 
much worse than a true planet fit. These configurations 
are therefore excluded. 

For blends involving a background/foreground star 
transited by a larger planet, there is in principle a wide 
range of allowed masses (spectral types) for the sec- 
ondary stars, as shown in Figure [TOl However, other con- 
straints available for Kepler- 18 strongly limit the number 
of these false positives. In particular, by comparing the 
predicted r — K s color of each blend against the measured 
color of the star fr om the KIC (r - K s = 1.723 ± 0.031; 
iBrownet al1[20lTI ). we find that many of the smaller- 
mass secondaries are ruled out because the blends would 
be much too red compared to the known color index of 
Kepler-18 (by more than 3c). Others are excluded be- 
cause the secondary star would be very bright (within one 
magnitude of the target, or in some cases even brighter 
than the target), and would have been noticed spectro- 
scopically, if unresolved in our AO or speckle imaging. 
This removes many but not all blends of this kind. 

For eclipsing binaries that are physically associated 
with the target (in a hierarchical triple star configura- 
tion) we find that the blend light curves invariably have 
the wrong shape to mimic a true transiting planet sig- 
nal, for any combination of stellar parameters for the sec- 
ondary and tertiary. Either the depth, duration, or steep- 
ness of the ingress/egress phases of the transits provide a 
poor match to the Kepler photometry. These scenarios 
are therefore all excluded. On the other hand, if we allow 
the tertiaries to be planets, then we do find a variety of 
secondary masses that can produce viable blends when 
transited by a planet of the appropriate size. The \ 2 map 
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TABLE 7 

Osculating Jaoobian Elements, at epoch 2455168.0 [BJD], for the 3-planet TTV Dynamical 

Solution. 



Period (days) To (days) ecosw esino; 



b 3.504674 ± 0.000054 266.276996 ± 0.005453 

c 7.641039 ± 0.000087 267.092502 ± 0.000262 0.000291 ± 0.000079 0.000173 ± 0.000233 

d 14.860509 ± 0.000148 269.174850 ± 0.000253 -0.000076 ± 0.000019 0.000516 ± 0.000450 



TABLE 8 

Masses and Densities of the Planets in the Kepler-18 System 



Method 






Kepler-18b 


Kepler- 18c 


Kepler-18d 








(K00137.03) 


(K00137.01) 


(K00137.02) 


MCMC (lightcurvc + RV) Solution 


(M a 


) 


13.4 ±5.8 


16.9 ±6.1 


29.9 ±8.8 


RV + transit time P and To 


(M a 


) 


12 ±5 


15 ±5 


28 ±7 


TTV dynamical model 


(M S 


) 


18 ±9 


17.3 ± 1.7 


15.8 ± 1.3 


TTV + RV dynamical model (adopted values) 


(M a 


) 


6.9 ±3.4 


17.3 ± 1.9 


16.4 ± 1.4 


density from adopted mass 


(g cm 




4.9 ±2.4 


0.59 ±0.07 


0.27 ±0.03 



for this general case appears in Figure ITTT and shows the 
range of radii permitted for the tertiaries, as well as the 
interval of secondary masses that yield suitable blends. 

The duration of a transit is set by the length of the 
chord traversed by the tertiary and the tangential veloc- 
ity of the tertiary during the event. The chord length, in 
turn, depends on the size of the secondary and the impact 
parameter. Therefore, the measured duration of a tran- 
sit provides a strong constraint on the allowed sizes for 
the secondary stars (or equivalently, their masses or spec- 
tral types). In the above BLENDER simulations we have 
placed no restriction on the orbital eccentricities of the 
star+star or star+planet pairs that can be blended with 
the target. When the orbits are permitted to be eccen- 
tric, the tangential velocity of the tertiary during transit 
can be significantly different from the circular case, and 
this allows a much larger range of secondary sizes than 



would otherwise be possible. In particular, chance align- 
ments with later-type stars transited by a planet near 
apoastron become viable as blends, and represent a good 
fraction of the false positives shown in Figures HU1 and [TP 
Some blends involving larger secondaries transited at pe- 
riastron can also provide acceptable fits. 

We note, however, that the period of K00137.03 is rela- 
tively short (3.5 days), and very large eccentricities, such 
as some of our simulated blends have, are unlikely as 
they would be ex pected to be damped by tidal forces 
(see [Mazch 2008j). Indeed, among binaries with main- 
sequence primary stars that are similar to Keplcr-18 
(solar- type, or later), all systems under 3.5 days have es- 
sentially circular orbi ts (see, e.g.. lHalbwachs et al.ll2003l : 
iRaghavan et al.ll2010[ ). We may take e = 0.1 as a conser- 
vative upper limit. Similarly, among the known transit- 
ing planets with periods of 3.5 days or shorter and host 
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Fig. 8. — Periodograms (fractional x 2 reduction as a function of 
frequency; [Zcchmcistcr & Kurstcr 200£|) of the O-C values shown 
in Figure [7] The dashed line shows the expected timescale from 
equation [JJ showing excellent agreement. The large tick marks 
show the Nyquist frequency l/(2P or t,it a i)> beyond which the peri- 
odogram holds no additional information. 
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Fig. 9.— Map of the Kepler-18b (K00137.03) \ 2 surface (good- 
ness of fit) for blends involving background eclipsing binaries with 
stellar tertiaries and arbitrary orbital eccentricities. The vertical 
axis represents the linear separation between the background bi- 
nary and the primary, cast for convenience in terms of the differ- 
ence in the distance modulus. Only blends within the solid white 
contour match the Kepler light curve within acceptable limits (3cr, 
where <r is the significance lev el of the x 2 difference compared to 
a transiting planet model; see[Frcssin ct al. 2011). Other concen- 
tric colored areas represent increasingly worse fits (4rr, 5rr, etc.), 
and correspond to blends we consider to be ruled out. Dashed 
green lines are labeled with the magnitude difference AKp be- 
tween the blended binary and the primary, and encompass the 
brightness range allowed by BLENDER (4 < AKp < 7). Blends with 
eclipsing binaries bright enough to be detected spectroscopically 
(AKp < 1 mag) are indicated with the hatched region below the 
solid green line, but are already ruled out by BLENDER. Similarly 
with the blue hatched areas that mark blends that are either too 
red or t oo blue compared to the measured color (see text and Fig- 
ure [lO}. When further constraining these blends to have realistic 
eccentricities (e < 0.1; see text), we find that all of them are ex- 
cluded by BLENDER. 

stars of any spectral type, none are found to have ec- 
centricities as large as e = 0.3. When constraining the 
false positives for K00137.03 to be within these eccentric- 
ity limits, we find that all background eclipsing binaries 
(stellar tertiaries) are easily excluded as they all require 
fairly eccentric orbits in order to match the observed du- 
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Fig. 10. — Similar to Figure[9](and with the same color scheme) 
for Kepler-18b (K00137.03) blends involving background or fore- 
ground stars transited by a larger planet. The blue hatched regions 
correspond to blends that are too red (left) or too blue (right) com- 
pared to the measured r — K s color of Kepler-18, and are thus ruled 
out. When blends are restricted to realistic orbital eccentricities 
(e < 0.3), many of the later- type secondaries are excluded (dashed 
3ct contour) The combination of the brightness (green hatched 
area) and color constraints leaves only a reduced area of parameter 
space ( "Allowed region" ) where blends are a suitable alternative to 
a transiting planet model. All of these scenarios have AKp < 7.0 
(dashed green line). 




0.4 0.6 O.f 

Secondary mass (M ) 

Fig. 11— Similar to Figure QJ] for Kepler-18b (K00137.03) for 
the case of hierarchical triple systems in which the secondary is 
transited by a planet. After taking into account the constraints on 
the r — K s color and brightness (blue and green hatched regions, 
respectively), only secondary stars with M2 < 0.5 Mq lead to blend 
light curves that match the observations. Further restriction to 
realistic orbital eccentricities (e < 0.3; see text) leads to a slightly 
smaller 3-cr contour (dashed). 

ration. Additionally, the numbers of blends involving 
star+planet pairs in the foreground/background or in hi- 
erarchical triple systems are considerably reduced when 
restricting the eccentricities, although many remain that 
we can not rule out. In the following we assess their fre- 
quency, and compare it with the expected frequency of 
transiting planets. 

7.1. Validating K00137.03 

The a priori frequency of stars in the background or 
foreground of the target that are orbited by a transiting 
planet and are capable of mimicking the K00137.03 sig- 
nal may be estimated from the number density of stars 
in the vicinity of Kepler-18, the area around the target 
in which such stars would go undetected in our high- 
resolution imaging, and the frequency of transiting plan- 
ets with the appropriate characteristics. To obtain the 
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number density (stars per square degree) we make us e 
of the Galactic structure models of iRobin et al.l (|2003f ). 
and we perform this calculation in half-magnitude bins, 
as shown in Table [TJ For each bin we further restrict 
the star counts using the constraints on the mass of the 
secondaries supplied by BLENDER (see Figure [TO]), and the 
eccentricity limit for transiting planets discussed above 
(e < 0.3). These mass ranges are listed in column 3, 
and the resulting densities appear in column 4. Bins 
with no entries correspond to brightness ranges excluded 
by BLENDER. The maximum angular separation (p maK ) at 
which stars of each brightness would escape detection in 
our AO/speckle imaging is shown in column 5 (see Ta- 
bleland Section l3"T2"j) . The result for the number of stars 
in each magnitude bin is given in column 6, in units of 
1CT 6 . 

To estimate the frequency of transiting planets that 
might be expected to orbit these stars ( and lead to a 
false po sitive) we rely on the results from iBorucki et all 
(|2011bf ). who reported a total of 1235 planet candidates 
among the 156,453 Kepler targets observed during the 
first four months of the Mission. These signals have not 
yet been confirmed to be caused by planets, and there- 
fore remain candidates until they can be thoroughly fol- 
lowed up. However, the rate of false positives in this 
sample is expected to be quite small (~10% or less; see 
iMorton fc Johnso"nll2011| ). so our results will not be sig- 
nificantly affected by the assumption that all of the can- 
didates are planets. W e further assume that the census of 
IBorucki et al.l ()2011bl ) is largely complete. After account- 
ing for the additional BLENDER constraint on the range of 
planet sizes for blends of this kind (tertiaries of 0.32- 
1.96 -Rjup), we find that the transiting planet frequency 
is /planet = 374/156,453 = 0.0024. Multiplying this fre- 
quency by the star counts in column 6 of Table Q] we ar- 
rive at the blend frequencies listed in column 7, which are 
added up in the "Totals" line of the table (12.7 x 10~ 6 ). 

Next we address the frequency of hierarchical triples, 
that is, physically associated companions to the tar- 
get that are orbited by a larger transiting planet able 
to mimic the signal. The rate of occurrence of this 
kind of false positive may be estimated by consider- 
ing th e overall frequency of binary stars (34% accord- 
ing to iRaghavan et al.l 120101 ) along with constraints on 
the mass range of such companions and how often they 
would be orbited by a transiting planet of the right size 
(0.43-0.53 i?j up ; see Figure ITT|) . The mass constraints 
include not only those coming directly from BLENDER, 
but also take into consideration the color and bright- 
ness limits mentioned earlier. We performed this cal- 
culation in a Monte Carlo fashion, drawing the sec- 
ondary stars from the m ass ratio distribution reported by 
IRaghavan et al.l (|2010[) . and the transiting planet eccen- 
tricities from the actual distribution of known transiting 
planets (http://exopla net.eu/[ ), with repetition. Planet 
frequencies in the appropriate radius range were taken 
as before from IBorucki et al.1 (|2011bf ) . The result is a 
frequency of hierarchical triples of 4.4 x 10 -6 , which we 
list at the bottom of Table [T] 

Combining this estimate with that of back- 
ground/foreground star+planet pairs described 
previously, we arrive at a total blend frequency of 
BF = (12.7 + 4.4) x 10" 6 w 1.7 x 10~ 5 , which represents 
the a priori likelihood of a false positive. From a 



Bayesian point of view analogous to that adopted to 
validate previous Kepler candidates, our confidence in 
the planetary nature of K00137.03 will depend on how 
this likelihood compares to the a priori likelihood of a 
true transiting planet, addressed below. 

The blend frequency of 1.7 x 10 -5 corresponds to false 
positive scenarios giving fits to the Kepler photometry 
that are within 3cr of the best planet fit. We use a sim- 
ilar criterion to estimate the a priori tra nsiting planet 
frequen cy by counting the KOIs in the IBorucki et al.l 
(|2011bl ) sample that have radii within 3a of the value 
determined from the best fit to K00137.03 (R p = 2.00 ± 
0.10 i?e). We find 284 that are within this range, giving 
a planet frequency PF = 284/156,453 = 1.8 x 10" 3 . 

This estimate does not account for the fact that the 
geometric transit probability of a planet is significantly 
increased by the presence of additional planets in the 
system (Kepler-18c and Kepler-18d in this case), given 
that mutual inclination angles in systems with multiple 
transiting plane ts have been found b e relatively small 
(typically 1-4°; ILTssauer et al.ll201l"bT) . Furthermore, a 
planet with the period of K00137.03 would be interior to 
the other two, further boosting the chances that it would 
transit. To incorporate this coplanarity effect, we have 
developed a Monte Carlo approach, described fully in 
Appendix El in which we simulate randomly distributed 
reference planes and inclination dispersions around this 
plane, from which a weighted distribution of the inclina- 
tion with respect to the line of sight for a third planet 
is calculated. Inclination angles relative to the random 
reference pl ane are assumed to fo llow a Rayleigh distri- 
bution (see iLissauer et aT1l2011bl) . Although the known 
planets carry some information on the inclination disper- 
sion, the probability of transit still depends somewhat 
on the allowed range of dispersion widths. When the 
assumed prior for the inclination dispers ion is uniform 
up to 4°, following Lis sauer et al.1 ()2011bl ). we find that 
the flatness of the system results in a very significant 
increase in the transit probability for K00137.03 from 
11.7% to 97%. To be conservative, we adopt a larger 
range of possible inclination dispersions from 0° to 10°, 
motivat ed by the upper lim it from the similar Kepler-9 
system (Hol man et al.ll2010l ). With this prior, the transit 
probability for K00137.03 becomes 84%, or an increase 
by a factor of ^7 over the case of a single transiting 
planet. 

Thus, the likelihood of a planet is more than 700 times 
greater (PF/BF = 0.013/1.7 x 10~ 5 « 700) than that 
of a false positive, which we consider sufficient to vali- 
date K00137.03 as a true planet with a high degree of 
confidence. We designate this planet Kepler- 18b. We 
note that our planet frequen cy calculation a ssumes the 
1235 candidates cataloged by IBorucki et al.l (|2011bl ) are 
all true planets. If we were to suppose conservatively 
that as many as 50% are false positives (an unlikely 
pro position that is also inconsistent with oth er evidence; 
see lBorucki et aLl^OllbUHoward et al.ll2010D . the planet 
likelihood would still be ~350 times greater than the like- 
lihood of a blend, implying a false alarm rate sufficiently 
small to validate the candidate. 

8. PHYSICAL PROPERTIES OF THE PLANETS 

The Kepler-18 system consists of two low-density 
Neptune-mass planets near a 2:1 mean motion resonance 
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TABLE 9 

Blend frequency estimate for K00137.03. 



Blends Involving Planetary Tertiaries 



Kp Range 


AKp 


Stellar 


Stellar Density 


pmax 


Stars 


Transiting Planets 


(mag) 


(mag) 


Mass Range 


(per sq. deg) 


(") 


(xlO- 6 ) 


0.32-1.96 fljup, /pi^et = 0.24% 






(AT©) 








(xlO" 6 ) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


13.5-14.0 


0.5 












14.0-14.5 


1.0 












14.5—15.0 


1.5 


0.87—1.40 


862 


0.08 


1.34 


0.003 


15.0-15.5 


2.0 


0.82-1.40 


1377 


0.11 


4.04 


0.010 


1 5 5-1 6 


2.5 


0.78-1.40 


2094 


0.13 


8.58 


0.021 


16.0-16.5 


3.0 


0.72-1.40 


3053 


0.16 


18.9 


0.045 


16.5-17.0 


3.5 


0.43-1.40 


4341 


0.20 


42.1 


0.101 


17.0-17.5 


4.0 


0.43-1.40 


5873 


0.24 


82.0 


0.197 


17.5-18.0 


4.5 


0.43-1.32 


7599 


0.32 


189 


0.454 


18.0-18.5 


5.0 


0.43-1.25 


9399 


0.40 


365 


0.876 


18.5-19.0 


5.5 


0.43-1.09 


10819 


0.56 


822 


1.973 


19.0-19.5 


6.0 


0.43-1.01 


11988 


0.64 


1190 


2.856 


19.5-20.0 


6.5 


0.43-0.92 


12585 


0.80 


1952 


4.685 


20.0-20.5 


7.0 


0.43-0.61 


3373 


0.88 


633 


1.519 


20.5-21.0 


7.5 












21.0-21.5 


8.0 












Totals 










5308 


12.7 



Blend frequency from hierarchical triples (see text) = 4.4 X 10 



Total frequency (BF) = (12.7 + 4.4) X 10~ 6 = 17.1 X 10~ 6 



Note. — Magnitude bins with no entries correspond to brightness ranges in which BLENDER excludes all blends. 



and an inner super-Earth-size planet. Its architecture 
bears a strong resemblance to Kepler-9, except that the 
Kepler- 18 system is less compact and its planets are less 
dense. The use of the observed transit times as well 
as the radial velocity data in the dynamical model of 
the Kepler-18 system places tight limits on the allowed 
planetary masses. We adopt these values as our best de- 
termination of the masses of the transiting planets in 
the Kepler-18 system. The last line in Table [1] gives 
the planet densities, computed from the final adopted 
planet masses. The TTV measurements together with 
the radial velocities restrict the masses of Kepler-18c 
(17.3±1.9M e ) and Kepler-18d (16.4±1.4M e ) to be sim- 
ilar to each other. Both are slightly lower than the mass 
of Neptune. Their radii however are 40% and 80% larger 
than Neptune respectively, giving them bulk densities of 
0.59±0.07g cm" 3 and 0.27±0.03g cm" 3 , which are only 
0.36 and 0.16 that of Neptune. The mass of Kepler-18b 
from the joint dynamical solution to the transit times 
and the RV measurements, is 6.9 ± 3.4M0. With its 
"super-Earth" size radius of 2.00 ± O.lOi?®, the density 
of this inner planet in the syste m is 4.9 ± 2.4g cm" 3 . 

Using the meth o ds des cribed in lMiller et al.l (|2009( ) and 
iMiller fc Fortnevl (|20ll , we have modeled the thermal 
evolution and interior structure of the two "Neptune- 
class" planets, Kepler-18c and Kepler-18d. Both planets 
have inflated radii compared to Uranus and Neptune, 
which points to two effects. The first is that the high 
incident flux slows their contraction. The second is that 
the mass fraction of heavy element within these two plan- 
ets is lower than that of Uranus and Neptune, meaning 
the mass fraction of H-He gas is larger. Uranus and Nep- 
tune are ~ 80-90% heavy elements (e.g ., water and rock) 
by mass (|Fortnev fc Nettehn ann 2010J), while below we 
show that the heavy element mass fractions of Kepler- 18c 



and Kepler-18d are somewhat lower than these values. 

Thermal evolution/contraction models are constrained 
such that the radius of each planet m ust be reproduced 
at the system's estimated age. As in IMiller fc Fortnevl 
(|2011l) . all relevant uncertainties are accounted for. 
These include uncertainties in the age of the system, the 
semimajor axes, masses, and radii of the planets, and the 
distribution of the heavy elements within each planet. 
We do not include an internal heating contribution due 
to eccentricity damping in either planet, as this power 
sourc e is expected to fall off as a -15 / 2 (Jackso n et al.l 
12008(1 and, furthermore, the eccentricities suggested by 
the TTV solutions are quite small. A 50-50 by mass 
ice-rock equation of state is used for the heavy elements. 
We find a heavy element mass of 13.5±1.8M® in Kepler- 
18c (~ 80% of the planet's mass) and 10.1 ± 1.4M® in 
Kepler-18d (~ 60% of the planet's mass). Kepler-18c is 
clearly more "metal-rich" than Kepler-18d. Both planets 
are consistent with a core-accretion formation scenario 
in which ~ IOMq of heavy elements gravitationally cap- 
tures an envelope of H-He gas. This envelope itself may 
be enhanced in heavy elements, as is inferred for Uranus 
and Neptune. 

Planet Kepler-18d, with a large radius of nearly 7i?©, 
may point to a population of Neptune-mass exoplanets 
having relatively low heavy element mass fractions and 
radii approaching that of the gas giant regime. They 
appear very similar to HAT-P-26b, a 4.2 day planet or- 
biting a cooler Kl-dwarf. The formation and evolution 
of such lower-densit y Neptune-class pla nets was recently 
studied in detail bv lRogers et al.l f|201 lh . They find that 
modestly more massive H-He envelopes than found on 
Uranus and Neptune (leading to larger planetary radii 
and lower densities) may be a common outcome of the 
core-accretion planet formation process. 
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Fig. 12. — The mass-radius diagram for transiting exoplanets (circles) and Solar System planets (stars). Left panel - Known transiting 
exoplanets (open circles) and Kepler-18b, c and d (filled circles ) , wit h curves of theoretical relations for cold pure hydrogen, water, 
rock (silicate), and iron (after Scagcr ct al. (2007); Fortncy ct al. (2007)). Two detailed models with mixtures and surface temperatures 
appropriate for Kepler-18b and Kepler-18c are also shown (IMillcr & Fortncy 2011). Right panel - Zoom on the smaller planets, with 
curves of constant mean density and curves of detailed interior models of constant composition. The theoretical models are (from top to 
bottom): "Neptune-clas s" models with hydrogen/helium envelopes and 50-50 ice-rock (by mass) cores from thermal evolution calculations 
(IMiller fc Fortncy 2011]) - the labeled H/He fractions are averages - if all metals are in a core, these fractions will be 31% and 16%, 
respectively; if mixed partially, they become 38% and 23%, respectively, for the curves shown; the '50% water' models have compositions 
of 44% silicate mantl e and 6% iron core, and the nominal 'Earth-like' composition with terrestrial iron/silicon ratio and no volatiles are by 
IValencia et~aH J200r3) and Zeng fc Sass elov (2011, subm.) The maximum mantle stripping limit (maximum iron fraction, minimum radius) 
was computed by IMarcus et ahl d2010l ). All these model curves of constant composition are for illustration purposes, as the degeneracy 
between compos ition mixture s and mean density is sign i ficant in this part of the m a ss-radius diagram. The data for the exopl anets 
were taken f rom IQueloz et al.l J2009) , Charbonncau ct al. (2009), Hartman ct al. (2011), Batalha ct al. (2011), Liss auer et al.l d2011al ) and 
IWinn et al.l fl2011D . We note that new analysis of CoRoT-7b lIHatzes et al.1120111) places it at a similar mass to Kepler-18b, and similar 
high-density composition to Kepler-lOb. The three unmarked exoplanets surrounding Kepler-18b in the diagram are (in increasing mass) 
CoRoT-7b, GJ 1214 b, and 55Cnce. 



The inner, 3.5-day period planet Kepler-18b, is a 
super-Earth that requires a dominant mixture of water 
ice and rock, and no hydrogen/helium envelope. While 
the latter cannot be excluded simply on the basis of 
the planet's mass and radius, the evaporation timescale 
for a primordial H/He envelope for a hot planet such 
as Kepler- 18b is much shorter than the old age de- 
rived for the Kepler-18 system, and such a H/He en- 
velope should not be present. Thus, despite its lower 
equilibrium temperature, Kepler-18b rese mbles 55 Cnc e 
and C oRoT-7b (a s originally measure d by IQueloz et al.1 
( 2009) though the Hat zes et al.l ([201 ll ) re-analysis makes 
CoRoT-7b very simila r to Kepler- 10b). Kepler-18b, to- 
gether with 55 Cnc e (jWinn et al.1 120111 ). are likely our 
best known cases yet of water planets with substantial 
steam atmospheres (given their high surface tempera- 
tures). 

It is interesting to compare the three transiting planets 
in Kepler-18 in terms of their apparent compositions and 
orbital sequence. Kepler-18 reinforces a pattern already 
seen in Kepler-11, and with less confidence in Kepler- 10 
and Kepler-9. Namely, inner planets are denser, though 



not always by very much - compare Kepler- 18b vs. c & 
d, and Kepler-llb,c vs. d, e, & f (Figure 12). It remains 
unclear at present whether this reflects a density gradient 
at formation or could be accomplished by evaporation 
later. 

Kepler was competitively selected as the tenth Discov- 
ery mission. Funding for the Kepler Mission is provided 
by NASA's Science Mission Directorate. We are deeply 
grateful for the very hard work of the entire Kepler team. 
This research is based in part on observations made with 
the Spitzer Space Telescope, which is operated by the 
Jet Propulsion Laboratory, California Institute of Tech- 
nology under a contract with NASA. Support for this 
work was provided by NASA through an award issued by 
JPL / Caltech. Some of the data presented herein were ob- 
tained at the W. M. Keck Observatory, which is operated 
as a scientific partnership among the California Institute 
of Technology, the University of California, and the Na- 
tional Aeronautics and Space Administration. The Keck 
Observatory was made possible by the generous financial 
support of the W. M. Keck Foundation. 
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APPENDIX 

INCORPORATING COPLANARITY IN MULTI-TRANSITING SYSTEMS TO ESTIMATE TRANSIT PROBABILITY 

In a system with N transiting planets, what is the geometric probability that an additional planet (called "planet 
iV+1") with a given period would transit, taking into account the fact that planetary orbits are expected (or observed) 
to be nearly coplanar? By incorporating the effect of coplanarity, the probability is significantly increased for additional 
planets to transit in known transiting planet systems compared to isotropically distributed planets. 

In a one-planet, one-candidate system, N = 1 and there is no t much that can be done unless a prior assumption 
for the true mutual inclinatio n is taken (Bcatty & Seager] 120101 ) . While in some cases, the true mutual inclination 
can be directly measured f see iRagozzine fe Holmanl ([2010) for various methods), generally it will have to be inferred 
from population statistics. iLissauer et al.l (|2011bl ) estimate that the inclination distribution of short-period planetary 
systems seen by Kepler ranges from 1-4°. This range could be used as a prior, both in the case of N = 1 and in the 
case of higher multiplicities. Additionally, with careful analysis, the absence of transit timing and duration variations 
can be used to p ut an upper limit on mutual inclinations, as in the case of Kepler-9 where the mutual inclination must 
be less than 10° ()Holman et al.ll2010[) . 

However, systems with two or more transiting planets, such as Kepler-18, have evidence of being thin without direct 
reference to the greater population of Kepler multiples. Furthermore, the inclinations of the known planets themselves 
give information on both the location of the reference plane of that system (i.e., the Laplace plane) and the typical 
inclination dispersion. If both planets have the same inclination slightly different from 90°, as in the case of Kepler-18c 
and d, then this suggests that the reference plane for this system is something like the average plane between the two 
planets and that the inclinations with respect to that plane a re likely quite small. In this situation, candidates with 
periods smaller than the inner planet are quite likely to transit (Ragozzine & Holman 2010). If both planets have quite 
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different inclinations, as in the case of Kepler-lOb and Kepler-lOc (Batalh a et al J 120 111) , this implies an uncertainty 
in the location of the reference plane and a non-flat system, and the probabilit y of additional planets transiting is 
enhanced over completely isotropic systems, but not as much (jFressin et al.ll201ll ). 

The goal is to quantify this effect in a natural way that uses all the available information. This is done using 
a Monte Carlo simulation that randomly generates reference planes, inclination dispersions, and nodal angles. The 
Monte Carlo trials that result in systems that match well the observed inclinations of the N planets are given higher 
weight than the vast majority of trials that do a very poor job. The line-of-sight inclination of the additional planet is 
also calculated in each Monte Carlo trial, allowing for the generation of a final weighted distribution. The assumptions 
made in this calculation are that exoplanet reference planes are, a priori, randomly oriented in an isotropic way and 
that the inclinations (with respect to the reference plane) of all the planets in a system are faithfully represented by 
a single Rayleigh distribution. The latter assumption does not account for the possible anticorrelation between size 
and inclination that can occur in dynamically thermalized systems, i.e., larger planets may be more well-aligned than 
smaller planets. With the Monte Carlo nature of the calculation and the error bars on inclinations that are typically 
produced by Kepler, the importance of a possible anticorrelation is reduced. 

In practical terms, the Monte Carlo simulation draws a random inclination (i r ) for the reference plane (uniform in 
cosi r ). The width of the Rayleigh distribution crj is drawn randomly from a distribution that is assumed a priori. 
For each Monte Carlo trial, a random inclination is drawn from the Rayleigh distribution as is a random nodal angle 
fo r each of the N + 1 planets; the nodal angles are chosen uniformly between and 2ir. Using the method in eq. 2 
of [Ra gozzine &: Holmanl (|2010h . the on-the-sky inclinations of the N planets are calculated. These are compared to 
the observed inclinations by computing the standard Gaussian z-score, i.e., by calculating the number of standard 
deviations away the trial value is from the observed value. (This could be modified if the inclinations and errors of the 
known planets are either not known or are known not to be Gaussian.) The assigned weight for each planet is equal 
to the area under the Gaussian that has more extreme scores than the trial value, i.e., with a z-score of z — '*' c ~'° t ' , 

the weight is w = 1.0 — erf(-^l), where iuc is the calculated trial inclination, i \, s is the observed inclination with error 

i err , and erf is the standard error function. If the trial value of the inclination is exactly the same as the observed 
inclination, the weight is 1; if it is many standard deviations off, then the weight is essentially 0. The total weight 
for a Monte Carlo trial is the product of these weights over all the known N planets based on their inclinations. As 
expected, the weight increases when the reference plane is taken near the average of the known planetary inclinations 
with a Rayleigh width similar to the standard deviation between the known inclinations. 

Each of the Monte Carlo trials also calculates an inclination for planet N + 1, based on the same reference inclination 
and &i. Using the weights derived from the observed inclinations of the N transiting planets, the weighted distribution 
of the inclination of planet N + 1 can be created. To calculate the probability that planet N + 1 would be found 
transiting, the sum of the weights for those simulations that would have produced a transiting planet is divided by the 
sum of the weights for the entire Monte Carlo simulation. 

Despite the information on the inclination dispersion from the known planets that is included in the weighting, the 
answer depends somewhat on the prior assumption for o~t. For example, when applied to Kepler-18 with Kepler-18b 
as planet N + 1, this technique predicts a transit probability of 47%, 84%, and 97% when the prior on <Ji is uniformly 
drawn between 0° and 90°, 0° and 10°, and 0° and 4°, respectively. In every case, the near coplanarity and low impact 
parameter of the two outer planets significantly increase the probability that Kepler-18b is transiting. Changing the 
prior on the inclination of the reference plane does not affect the result. This method provides a quantitative way 
to estimate the increase in transit probability for planets in multi-transiting systems. When applied in combination 
with BLENDER and other techniques, it will allow multi-transiting systems to be validated more easily than singly- 
transiting systems. Note that the method described here only estimates the improvement in the probability of the 
planet hypothesis due to geometric const raints, and does not i nclude the additional effect that planets tend to be 
found in multiple systems, as discussed in lLissauer et al.l ()2011b( ). 



